UTTG-09-06 



A No- Truncation Approach to Cosmic Microwave Background 

Anisotropics 



Steven Weinberg* 
Theory Group, Department of Physics, University of Texas 
Austin, TX, 78712 

Abstract 

We offer a method of calculating the source term in the line-of-sight integral for cosmic 
microwave background anisotropies without using a truncated partial-wave expansion in the 
Boltzmann hierarchy. 



Electronic address: weinberg@physics.utexas.edu 



1 



I. Introduction 

Originally the Boltzmann equation for the photon distribution in cosmology [1] was solved 
numerically by expanding the components of the photon density matrix in a series of Legen- 
dre polynomials [2], but to get results that could be compared with observation of the cosmic 
microwave background this method requires the inclusion of hundreds or even thousands of 
partial waves, requiring hours or even days of computer time for each theoretical model. A 
great improvement was introduced with the suggestion to use instead a formal solution of 
the Boltzmann equation, in the form of a "line of sight" integral [3]. But this is still only a 
formal solution, in the sense that we still need to calculate source terms appearing in the 
integrand. These terms involve partial waves for the photon distribution up to £ = 2 for 
the scalar modes and £ = 4 for the tensor modes, and these of course are coupled to higher 
partial waves. In the original proposal of the "line of sight" method, and in the computer 
programs CMBfast and CAMB based on this method, these source terms are calculated 
numerically, by first finding an approximate solution of the Boltzmann equations for partial 
wave amplitudes. An accurate solution for the partial waves appearing in the source terms 
can be found by truncating the partial wave expansion at a sufficiently high value of £. In the 
latest version of CMBfast, the integrand for scalar modes is calculated using partial waves up 
to £ = 12, in which case one has to solve at least 26 coupled ordinary differential equations 
for the evolution of the partial wave amplitudes, not counting the equations needed to follow 
the evolution of the baryonic plasma, cold dark matter, neutrinos, and gravitational field 
components. For tensor modes, the source terms are calculated by solving the 22 differential 
equations for partial waves with up to £ = 10 for photons. The results for £ < 2 or £ < 4 
are used in this method to calculate the integrand of the line-of-sight integral, which then 
is used to calculate all the higher partial wave amplitudes measured in observations of the 
cosmic microwave background, up to values of £ over 1000. 

This article will present an alternative approach, which does not use partial wave ex- 
pansions to calculate the source terms, and hence obviates the need for any truncation of 
this expansion. Instead of a large number of coupled differential equations for the partial 

2 



waves, we have a single integral equation for the tensor modes, and a trio of coupled integral 
equations for the scalar modes (including one for the plasma velocity). Of course, integral 
equations are generally harder to solve numerically than differential equations (no routine 
for solving them is supplied by Mathematica), but in the case at hand the integral equations 
can be solved numerically by simple iteration. In this method, the calculation itself provides 
an immediate way of judging its own reliability — if the nth iteration agrees with the n — 1th 
iteration to a satisfactory degree of accuracy, one has a solution. In a sample calculation of 
the source term for the tensor modes, the results converge rapidly in just a few iterations. 

This paper concentrates in the next section on the calculation of the photon distribution, 
but a truncated partial wave expansion is also unnecessary for neutrinos. Indeed, it is 
already known [4] that the momentum distribution of massless neutrinos for a given metric 
perturbation can be calculated in terms of a simple line-of-sight integral, with no need to solve 
integral equations. CMBfast does not use this line-of-sight method for neutrinos, presumably 
because no one is interested in very high partial waves in the neutrino distribution, but to 
get good accuracy for the neutrino contribution to the energy-momentum tensor it carries 
the partial wave expansion to £ max = 25. In the last section of this paper the approach 
of reference [4], which dispenses with partial wave expansions, is extended to the case of 
massive as well as massless neutrinos, and to scalar as well as tensor modes. 

II. Photons 

First, some reminders about the Boltzmann equation for photons in cosmology. We will 
adopt a coordinate system in which the metric takes the form 

g>oo(x, t) = -1 , 0oi(x, t) = , ^(x, t) = a 2 (t) (5ij + hijixi, tfj , (1) 

where is a first-order perturbation. Weakly perturbed metrics will automatically be of 
this form in tensor modes, and can be put in this form for scalar modes by adopting a 
synchronous gauge. 

For our purposes, it is important to write the Boltzmann equation for the photon dis- 
tribution in a matrix form, rather than in the partial wave formalism in which it is usually 
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presented. The photon distribution is described by a polarization density matrix n J:, (x, p, t), 
defined so that if we measure whether photons have polarization in a direction e l rather 
than in an orthogonal direction, then the number of photons with polarization e l in a vol- 
ume rij d-Pi dx l of phase space at time t will be found to be gikgjiG k e l n l:l (x, p, t) Y\ m dp m dx m , 
with piii^ = 0. (The polarization of a photon with 3- momentum p { is described by a polar- 
ization vector e\ satisfying pit 1 = and g^e % e d = 1.) For small perturbations, this matrix 
can be put in the form 



n t3 (x.,p,t) = -n 7 ya(t)^g kl (^t)p k pi 

+ 5n ij (x,p,t) . 



g ij (*,t) 



£ r(x,t)g~ ? '(x,t)p fc p; 
g H (x,t)p k pi 



Here n 7 (p) is the equilibrium phase space number density 
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n 7 (p) 



exp [p/ksa{t)T{t 



(2) 



(3) 



(which is a time-independent function of its argument because in the era of interest T(t) oc 
l/a(t)), and 5n 1J is a small perturbation. This perturbation satisfies a linearized Boltzmann 
equation: 

dSn ij {x,p,t) p k d5n ij (x,p,t) 2d(t) 



+ 



+ 



dt ' a(t) dx k a(t) 



Sn l3 (x,p,t) 



Aa 2 (t) 
—u c (t) 5n lJ (x, p,t) + 



8tt J 



d% 



x [<5n u (x, ppi , t) - pipk 5n kj (x, pp x , t) - Pjp k Sn lk (x, pp x , t) 
+ PiPjPkPi 5n k \y.,ppi,t) 

UJ c (t) 



2a 2 (t) 



p k 5u k (x,t)n' (p) Sij-pipj 



(4) 



where p = ^fp^pli pk = Pk/p, 8u k is the streaming velocity of the baryonic plasma, and u c is 
the frequency with which a photon collides with electrons in the plasma. Instead of 5n % 3 , it 
is sufficient to consider the intensity matrix perturbation J^(x,p, t), defined by 



a 4 (t) pj (t) Jij (x, p,t) = a 2 (t) / 5n tj (x, pp, t) Anp 3 dp , 

J 



(5) 



where p 7 (£) = a~ 4 (t) / 47rp 3 n 7 (p) dp is the mean photon energy density. (This is all we need 
to calculate the photon contributions to the perturbations in the energy- momentum tensor.) 
To derive the Boltzmann equation for Jjj(x,p, t) we multiply Eq. (4) with 4np 3 and integrate 
over p = ^PiPi, and find 

dJij(x.,p,t) p k dJjj(yL,p,t) 
dt a(t) dx k 

+ p k Pihki(x, t) (Sij - pipj) 
= -u c (t) Jij(x,p,t) + / d 2 Pl 

x 



8tt 

Jij(x,pi,i) - PipkJkj{*,Pi,t) - PjpkJik(*,Pi,t) 



+ PiPjPkPi J k i(*,Pi,t) 



Hj PiPj 



p k Su k (x,t) . 



(6) 



We can also calculate the equation of motion of the plasma, using the perturbed momentum 
conservation equation: 

p B (t) d 



n 4 f d^p 

a(t)8u k (x,t)\ = --w c (t)p 7 (t)6uk(-x.,t) +w c (t)p 7 (t) J — J u (x,p,t)p k , (7) 



a(t) dt 

where ps is the unperturbed density of the baryonic plasma. These partial differential 
equations can be converted to ordinary differential equations by writing the perturbed metric 
as a Fourier integral 

h tJ (x,t) = J (Fpe^hifat) . (8) 
and looking for solutions in the form 



Jij(x,p,t) = J Ae iq 'Vy(q,p^) , 6ui(x.,t) = J d 3 q e iq-x <^(q,£) . 



(9) 



Then Eqs. (6) and (7) become ordinary differential equations, with d/dx l replaced with 
They have a formal solution in the form of a "line of sight" integral [3], which in the general 
case may be written 



rt / rt dt" rt N 

Jij(q,p, t) = ^ dt' exp I -iq • p - dt" u c (t") 



- PkPi (Sij - pipj) fe«(q, 



x 

+ 3uJc ^ (Jij(q, t') - PiPkJ kj (q, if) - pjPkJikiq, t') + PiPjPkPiJki(<l, If)) 



+ 2u c (t') [S i:j - pipj] p k Su k (q, If) 



+ JiM,pM) , 



(10) 



Aa(t) Jti 
+ <Jtti(q,*i) . 



dt 



„ co c (t") \ tu c (t')a(t>) 
R(t") ) R(t') 



£(q,0 



(11) 



We have here introduced the convenient abbreviations 

— Jij(q,p,t) , 

/(pp 
— Jkk(<l,P,t)Pi , 

and, as usual, 



(12) 
(13) 



R(t) 



3p B (t) 



4p 7 (t) • 

If we take the initial time t\ early enough so that photons are in local thermal equilibrium 
with the plasma at that time, then 

~5T(clM) 



T(h) 



+ p k Su k (q,t 1 ) 



(14) 



It is the calculation of the source terms Jij(q, t) and 2j(q, t) that concerns us in this 
paper. For this purpose, we now need to distinguish between tensor and scalar modes. We 
first consider tensor modes, which are computationally simpler. 

Tensor Modes 

For tensor modes the metric perturbation takes the form 



A=±2 



(15) 
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where /3(q, A) is a stochastic parameter; ey(g, A) is a polarization tensor for helicity A, with 
cji eij(q, A) = and e kk (q, A) = 0; and T> q (t) is the solution of the wave equation 

2 

V q (t) + Z-V q {t) + \v q (t) = 16nG Ti T q {t) , (16) 
a 



that does not decay while outside the horizon. Here ir q (t) is the coefficient of J2\ &ij(q, A)/3(q, A) 
in the Fourier transform of the tensor part of the anisotropic inertia tensor. (We are consid- 
ering times sufficiently late so that the decaying solution makes a negligible contribution to 
Sgij.) The quantity Jij(q, t) will then take the form of a corresponding sum over graviton 
helicities 

J ij {q,p,t)= /3(q,A) Jij(q,p,*,A) , (17) 

A=±2 

with Jij(q,t, A) ordinary c-number functions, not stochastic fields, satisfying the Boltzmann 
equation 

+ PkPi e k i(q, A) V q (t) fa - pipj) 
= -uj c (t) Jij(q,p,t,X) + 



2 



x 



Jij{q, t, A) - pip k J kj (q, t, A) - pjp k J ik {q, t, A) 



+ PiPjPkPiJki(q,t,X) , (18) 



where 



j2 a 

^■(q,t,A) = | Jij{q.,p,t,\) . (19) 

(The velocity perturbation <5-Uj is absent in tensor modes.) Furthermore, because Jy(q, t, A) 
for a given helicity A must be a linear combination of the polarization tensor components 
e k i(q, A) with the same A, while q k e k i(q, A) and e kk (q, A) both vanish, the only possible form 
of J7y(q, t, A) allowed by rotational invariance is just e^g, A) times some function of q = |q| 
and t This relation is conventionally written 

J ij (<l,t,\) = --e ij {q,\)*(q,t) . (20) 



To make contact with the notation used in the usual calculation of the source term \I/(g, t), 
we note that the intensity matrix perturbation may be written in the form 

Jij(q,p,t,X) = -(5ij - PiPj) PkPi e M (q, \)(aP (q,p ■ q,t) + A ( p\q,p ■ q, tj) 
+ (eij(q, A) -PiPke-kMA) -PjPkeik{q,\)+PiP j PkPie k i{q,\)) A ( p (q,p ■ q,t) . 

(21) 

(Here the superscript T stands for 'tensor,' while the subscript T stands for 'temperature.' 
The coefficients are chosen so that Ju is proportional to At, and the polarization is pro- 
portional to Ap.) A third term proportional to (g« — pi(p • q))(qj — pj(p • q))pkPi&ki would 
be allowed by symmetry principles, but is not generated by the Boltzmann equation. Using 
Eq. (21) in Eq. (18) yields separate Boltzmann equations for AP and A^P: 

-AP(q,p-q,t) + ta~\t)q-pA^(q,p-q,t) = -2V q (t)-uj^ , 

(22) 

^ A { P (q, p-q,t)+ia-\t)<l-p A { P (q,p-q,t) = -u c (t) AP (q, p ■ q, t) - u c (t) (q, t) . (23) 

To calculate the source term ^(q, t) in terms of partial waves, one first integrates Eq. (21) 
over p, and finds 

+ ((p-q) 2 + l(l-(p-q) 2 ) 2 ) A ( p(q,p.q,t) . (24) 
The functions AP and aP may be expanded in Legendre polynomials 

oo 

AP(q,p.q,t) = J2 i ~ e (^+ 1 ) P ^-P) A T r hq,t) (25) 
e=o 

oo 

AP(q,p.q,t) = Y.^W+l)PM-p)AP{q,t), (26) 

and then Eq. (24) reads [5]: 

*(?,*) = ^A^(9,0 + \Ag\q,t) + ^A ( P(q,t) - ?Ag2(9,0 

+ ^AP(q,t)--^AP(q,t). (27) 



As already mentioned, experience shows that to accurately calculate the partial wave am- 
plitudes up to £ = 4, which appear in Eq. (27), one needs to solve the Boltzmann equations 
for the partial wave amplitudes up to larger values of £, up to £ = 10. Once ^ is calculated 
in this way, we can calculate Jij(q,p, t, A) for very much higher values of £ by using the "line 
of sight" integral (10), which for the tensor modes gives 



Jij 



rt t rt df" rt \ 

q, p, t, A) = ^ dt' exp (-iq • p ^ - ^ dt" u c {t") j 

- PkPi (Sij - Pipj) e kl (q, A) t> q (t) 
-u> c (t') V(q,t')(eij(q, A) - PiP k e k j(q, A) - p j p k e ik {.qA) + PiPjPkPiekMA 



(28) 



(In tensor modes there is no perturbation to either the temperature or the velocity of the 
baryonic plasma, so Eq. (14) gives the initial value Jy(q, p, t\) = for an initial time t\ 
taken sufficiently early so that photons are in local thermal equilibrium with the plasma.) 

Here we suggest the alternative, of deriving an integral equation for ^(q, t) by simply 
analytically integrating Eq. (28) over p. Equating the coefficients of on both sides gives 
the integral equation 



V(q,t) 



ti 



dt' exp 



UJ c (t") dt" 



x 



tf(g,f) 



(29) 



Here K(v) and F(v) are the functions 



K(v) = j 2 (v)/v 2 , F(v) = j (v) - 2 3l (v)/v + 2 j2 (v)/v 2 



(30) 



Integral equations of this sort are harder to solve numerically than differential equations, 
because in a step-by-step calculation it is necessary to keep track of ^(t') for all t' < t in 
order to calculate ^(t). On the other hand, such equations can also be solved by simple 
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iteration, in which the numerical work is reduced to doing a few integrals. We calculate the 
nth approximation to \I/ by using the previous approximation vl/t™™ 1 ) in the second term 
in square brackets in Eq. (29), and start this calculation by taking the lowest approximation 
as just Eq. (29) with the second term in square brackets dropped: 

¥°\q,t) = -3jfV exp [- £ u c (t") dt"~\ V q (t')K ■ (31) 

Using \l/( ) for \l/ in Eq. (28) would give precisely the same photon distribution function that 
we would find if we assumed that photons remained unpolarized until the last scattering; 
subsequent iterations take account of the polarization produced in multiple scattering. When 
after n iterations we find that \&( n ) = \I/( n_1 ) to an adequate degree of accuracy, we have 
a solution of Eq. (29). This is in contrast with the usual truncation method, in which the 
accuracy of the calculation for a given maximum multipole order can only be judged by 
comparing with the results for a higher maximum multipole order. 

To test the convergence of this iteration procedure, we adopt the usual ACDM model, 
with cosmological parameters 

tt B h 2 = 0.0223 , tt M h 2 = 0.1262 , h = 0.732 . 

We calculate the collision rate u c by solving the kinetic equations for hydrogen recombination, 
with an un-ionized helium fraction Y = 0.26. The numerical calculations are done with q 
chosen so that the physical wave number q/a comes within the horizon at just the time of 
matter-radiation equality. The gravitational wave amplitude T> q (t) is calculated ignoring the 
effect of anisotropic inertia. Over the interesting time interval when the matter/radiation 
density ratio y increases from 2 to 4, which includes the time of recombination, the first 
iteration ^^(t) has roughly the same time- dependence as ^ q °\t), but is about 13% to 33 
% larger. On the other hand, after five iterations we get a result ^!^(t) that at worst (at 
y ~ 3) is within 0.3% of the previous iteration \I>( 4 ), and is much closer at other values of y. 
As a further illustration of the convergence of the iteration procedure, Figure 1 shows the 
first few iterations for several wave numbers and a broader range of values of y. 
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We will not go on here to use this source function in Eq. (28) to calculate the intensity 
matrix perturbation, because we are not proposing anything new in that part of the calcu- 
lation of the microwave background anisotropies, but only in the calculation of the source 
terms appearing in the integrand of the line of sight integral. The sample calculation de- 
scribed here shows that this is a practical method of calculating the source terms, as well as 
having the conceptual advantage of avoiding a more-or-less arbitrary truncation of a partial 
wave expansion. 

Scalar Modes 

For scalar modes, the metric perturbation in synchronous gauge takes the form 

hij(q, t) = 5ijA(q, t) - q^ B(q, t) . (32) 

Assuming the perturbation to be dominated by a single mode (presumably the adiabatic 
mode that does not decay while outside the horizon), the dependence of A(q, t) and B(q,t) 
on the direction of q is entirely contained in a stochastic factor f3(q) for this mode: 

A(q, t) = /%) A q (t) B(q, t) = /%) B q (t) . (33) 

In this case, the source terms (12) and (13) appearing in the integrand of the line-of-sight 
integral must take the form 



JuM = /%) 



(34) 



2i(q,t) =iP(q L )q i n(q,t) (35) 
Then, using Eq. (11) to evaluate the plasma velocity, the line-of-sight formula (10) becomes 

J^(q,p, t) = /%) dt' exp f -iq ■ p ^ — - - ^ dt" u c (t")j 



x 



(^■-M-) (A q (t')-(p- q fB q (t')) 

${q, t') (Si, - pipj) + -U(q, t') - pi(p ■ g)) (^- - pj(p ■ g)) j 
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+ - m) <p ■ » £ *• ^P«(^. n «p (- /I * 

+2w c (t') (<J - - p^) p k 5u k (q, *i) + 2(<% - ftp,- 



,cu c (t" 



R(t" 

T^TT-x +PkSu k {q,t 1 ) 
1 W 



(36) 



Conventionally, the intensity matrix perturbation for scalar modes is written 
Jij(q,P,t) = /%){^(A T ?) (g,g-p,t) - A { p\q, q ■ p, t)) - 



+ A P S) (g,g-p,t 



(ft - (9 • (& - (? • P)Pj) 
1 - (p • g) 2 



with A^ and Ap ; satisfying the Boltzmann equations 



jA). (q. f ,.l) = -w c (i)A£'(g,M) + 4^c(0 (1 -^)n(g,t) , 



(37) 



(38) 



A { T S \q,ti,t) + il-^r ) A ( T S \q,fi,t) = -u c (t)A { T S \q,^t)-2A(q,t) + 2q 2 ^B(q,t) 



At), 
3 

+ 3iu c (t) $(g, i) + -u c (t)(l - /i 2 )n(g, t) + Auj c (t)pi5 Ui ; . 

(S) (S) 

By expanding A T and A P in partial wave amplitudes 

00 

A?W,0 = £i^(2€+l)P,0x)A$(«z,f) 

00 

a^W,;) = Er^ + ijp^A^V), 



(39) 



(40) 
(41) 



and integrating Eq. (37) over the directions of p, one obtains well known expressions for the 
source functions in terms of the partial wave amplitudes with £ < 2: 

1 



$ = _ 

6 



TO ^PO AA T2 ^P2 



n 

ft 



A (S) , A (5) a(5) 



i P0 -r ^v T2 



4 P2 



i T1 



(42) 
(43) 
(44) 
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Instead, we suggest the alternative of deriving coupled integral equations for <3>, II, and 
CI by analytically integrating over p in Eq. (36). This gives 



where 



t) = J* dt' exp (- J* dt" w c (*")) 



2 {"•(*V,(« J f^) + ;n Wl (, J f*)} 

3 „ / f' <i(" \ C'',,, w 1 .((")o.((")„, ,„ / r*' ui c (t'")dt 



n(g, t) = jf* eft' exp jf* di" w c (t") 



x 



*M0{«< 9 , OA U SI + n (9 , Oft (, ^ ) } 



(it" 



^ dt" \ ft' uj c (t")a(t") 



i(t") J Jh 



R(t") 



' u c {t'") dt" 
R(t"') 



n(q, t) = jf dt' exp jf* dt" uj c {t") 



x 



3a; c (t / ) 
2 

3 ^ / ft dt" \ ft' , „ uj c (t")a(t") ... 



n> e (t"') g 
' i?(t w ) 



(45) 



(46) 



(47) 



Fx(t;) 



ji(v)/v - j (v) , 

ji(v)/v-j 2 (v) - j 2 (v)/v 2 + j 3 (v)/v(v) , 
ji(v)/v + j 2 (v)/v 2 - j 3 (v)/v , 



(48) 
(49) 
(50) 



13 



F 4 (v) = ji(v)-j 2 (v)/v, (51) 

F 5 (v) = -2j 2 (v)/v 2 + 5j 3 (v)-j 4 (v), (52) 

F 6 (u) ee j Q (v) - 2 3l (y)/v + 2j 2 (u) + 2j 2 (v)/v 2 - 5j 3 (v)/v + j 4 (v) , 

(53) 

F 7 (v) ee -2j 2 (^)/ V + j, 3 (^) , (54) 

F 8 (u) ee 3j 2 (v)/v-j 3 (v) , (55) 

F 9 (u) ee -j 1 (v) + 3j 2 (v)/v-j 3 (v) , (56) 

F 10 (u) ee 2j 1 ( U )/u - 2j 2 (v) , (57) 



(We have dropped the terms involving the plasma temperature and velocity perturbations 
at the initial time t±, because as long as the collision rate is much larger than the expansion 
rate the photon distribution is given in terms of ST/T and 8u by an equilibrium formula like 
Eq. (14), and as long as the perturbation is outside the horizon 5T/T and 8u grow in the 
adiabatic mode, so that if ti is taken sufficiently early in the era when the collision rate is 
large and the perturbation is outside the horizon, the initial-data terms become negligible.) 
One may again expect that these coupled integral equations may be solved by iteration, as 
done for the tensor modes. 

III. Neutrinos 

In the absence of collisions, the Boltzmann equation for either massless or massive neu- 
trinos for the metric (1) reads 

dn p l dn dn jpp k dgjt . . 

~dt + ]£Ihf = ~d^~^~ ' 

Here n(x, p, t) is the phase space density of neutrinos, regarded as a function of the compo- 
nents x\ Pi and the time t, while p l and p° are functions of the x l and t as well as the p iy 
given by p % = g tj pj and p° = sJg^PiPj + m 2 . (The derivation of Eq. (58) is given in Ref. [4] 
for massless neutrinos, but precisely the same derivation applies also for massive neutrinos.) 
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For weak perturbations, we write 

ra(x, p, t) = n v (ait) yj 'gv (x, t)p(p^ + <Jn(x, p, t) , (59) 
where n„ is a time-independent function of its argument 



n„(p) 



cxp 



0? 2 + a 2 (ti)m 2 \ 
k B a(t)T(t) 1 ' ' 



(60) 



(27T) 

Here 5n is a small correction, representing the dynamical rather than purely geometric 
effect of metric perturbations on the neutrinos, and t 1 is the time the neutrinos went out of 
equilibrium with the baryonic plasma. (Probably all neutrinos have masses much less than 
ksTiti), in which case the term a 2 (ti)m 2 in the square root can be neglected.) Expanding 
to first order in 5n and 8gij = a 2 hij, after many cancellations this gives 

dSn + d5n Pi _ n '^ PkPk \ KiVKVi (gi) 

dt dx i a^pkPk + aFrrP 2^/p k p k kl k 1 ' 

We again write hij as a Fourier integral (8), and seek a solution in the form 

5n(x, p,t) = J d 3 q e iq - x 5n(q, p, t) . (62) 

Then Eq. (61) becomes an ordinary differential equation, with d/dx l replaced with iq^. 
Instead of solving this equation numerically with a truncated partial wave expansion, it can 
be solved analytically, as a line of sight integral 

dn(q, p, t) = / dt exp — iq ■ p / 

2^/pm Jti \ h> a(t")^/p kPk + a 2 (t")m' 2 

x h k i(q,t')p k pi . (63) 

where t\ is any time taken early enough so that 5n(x, p, t\) is negligible. 

For the foreseeable future, the neutrino distribution will be important only in calculating 
the neutrino contribution to the energy-momentum tensor. The first-order perturbation to 
the contribution of each species of neutrino or antineutrino to the mixed components of the 
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energy-momentum tensor takes the simple form 



*Z*-(x,f) = 



1 



a\t) 
1 



Yl dp k £ra„(x,p,t) 



PiPj 



\jvWk + a?{t)m 2 
\\ dp k \ Sn v (x,p,t)pj , 



\k=l 
< 3 



vfe=l 



^°o(x, f) = / (il d Pk) ^( x > P> \jvkPk + a 2 {t)m 2 



(64) 

(65) 
(66) 



Kt) - v, , 

with all other first-order contributions nicely cancelling. In the tensor mode, for which 
qihij = ha = 0, the only non-vanishing perturbations are to the space-space components 
of the energy-momentum tensor, which are needed in calculating the viscous damping of 
gravitational waves [4]. These take the form 



d 3 qe^ x j Q 4np 5 n' u (p)dp 
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x f dt' K \qp C 

JU \ Jt> n (+"' 



dt" 



hij{q, if) 



t' a{t")Jp 2 + a 2 {t")m 2 / Jp 2 + a 2 (t')m 



(67) 



where K(v) = j2(v)/v 2 . For scalar modes, we get contributions to all components, of the 
form 

4 

5T l vj = 5ij5p u + didj-Kv , 5T®j = -p u dj5u u , 5T° = -5p u , (68) 

in which Sp u , ir u , 5u v and 5p v may be identified as the pressure perturbation, scalar anisotropic 
inertia, velocity potential, and energy density perturbation, respectively, for a given neutrino 
species. Inserting Eqs. (32) and (33) in Eq. (63) and then using the result in Eqs. (64)-(66) 
gives 

K ' M(t)J ^Jo ^f + a 2 {t)m 2 

dt" 



x 



fdt' A q (t')F n [qp f 

-q 2 B q (t')F 12 [qp f 

\ Jf n +" 



)^/p 2 + a 2 (t")m 2 i 
dt" 



*' a(t")Jp 2 + a 2 (t")m 2 



(69) 
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p 2 + a 2 (t)m 2 
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-B q (t')F 5 {qp ( 

Jt ' a(t")Jp 2 + a 2 {t")m 2 



A I f roo 

-p u (t)5u v (x, t) = -— . y d 3 g e^ x /3(q) y 47rp%'» dp 

X /'dt' q^AJt'folqp [* , dt 

Jti V *' a(t")^/p 2 + a 2 (t")m 2 / 
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* a(t") x /p 2 + a 2 {t")m 2 
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a(t")^Jp 2 + a 2 (t")m 2 



where 



Fn(v) =ji(v)/v , F 12 (v) = j 2 {v)/v 2 - j 3 (v)/v , F 13 (v) = ji(v)/v - j 2 (v) 



(70) 



(71) 



(72) 



(73) 



It is only in the case m = that the integrals over neutrino energies can be done separately 
from the integrals over time, and give results simply proportional to p v . 

The contribution of photons to the perturbations in the energy momentum tensor can be 
calculated in a similar way, by using Eqs. (64)-(66) with a 2 Sn u in place of Sn (and of course 
with m — 0), taking the integral of 5n tl over photon energies from the line of sight integrals 
(28) or(36) for tensor or scalar modes, respectively. 

Added Note: After the preprint of this paper was first circulated, I learned that similar 
suggestions regarding the calculation of the photon source terms have been made in the 
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preprint of a paper by D. Beskaran, L. P. Grischchuk, and R. G. Polnarev, gr-qc/0605100 



Their equation (61) is the same integral equation for the source terms of tensor perturbations 
to the cosmic microwave background as presented here in equation (29). However, for scalar 
modes they do not give an integral equation for the baryonic plasma streaming velocity, and 
so instead of the three coupled integral equations found here, they give two coupled integral 
equations, in which the plasma velocity as well as the gravitational field perturbations appear 
as inputs. 

I am grateful to Eiichiro Komatsu for frequent helpful conversations, and to Raphael 
Flauger for preparing the figure and pointing out some typographical errors. This material 
is based upon work supported by the National Science Foundation under Grant Nos. PHY- 
0071512 and PHY-0455649 and with support from The Robert A. Welch Foundation, Grant 
No. F-0014. 
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Figure 1: Iterative Solution of Equation (29) for the Source Function \E'. Short dashes 
indicate the zeroth iteration (31); longer dashes indicate the first iteration; and further 
iterations are indistinguishable from the solid curve, which therefore represents the solution. 
Here the source function is calculated for a gravitational perturbation V q {t) whose value V° 
before horizon entry is unity; for any other initial condition, \1/ should be multiplied by the 
value of V°. The quantity n is the wave number q in units of the wave number that just 
comes into the horizon at matter-radiation equality, and y is the Robertson-Walker scale 
factor, in units of the scale factor at matter-radiation equality. The calculation was done 
by R. Flauger using the electron number density calculated with the program Recfast [S. 
Seager, D. D. Sasselov, and D. Scott, Astrophys. J. 523, LI (1999); Astrophys. J. Suppl. 
128, 407 (2000)] and taking Q M h 2 = 0.133, Vt B h 2 = 0.02238, h = 0.72, T = 2.725K, and 
Y He = 0.24. 
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